差异分析及功能注释（上）

课程笔记

粉丝：有单细胞线上课程吗？

小编：什么？ 我们的单细胞转录组分析线上课程已经上线好久了，你们竟然都不知道吗，每篇推文后面的课程推荐没人看的吗，小编已哭晕在厕所

好了，戏演完了，下面郑重介绍下我们的单细胞线上课程：（详情戳下方链接） 

这个课程笔记栏目记录了学员们学习单细胞转录组课程的学习笔记

希望大家能有所收获！

作者 | 单细胞天地小编 刘小泽

课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
这次会介绍如何进行差异分析及功能注释，以及比较不同分组。对应视频第三单元9-11讲

前言

前面得到的6个发育时期和4个分群，而且还可视化了一些marker基因，那么现在就要对这4群细胞进行差异分析

将对应文章这张图：

1 使用monocle V2进行差异分析

1.1 准备表达矩阵和分群信息

 1rm(list = ls()) 
 2options(warn=-1) 
 3options(stringsAsFactors = F)
 4source("../analysis_functions.R")
 5
 6 # 差异分析一般都要求使用count矩阵
 7load('../female_count.Rdata')
 8
 9# 6个发育时期获取
10head(colnames(female_count))
11female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1)
12names(female_stages) <- colnames(female_count)
13table(female_stages)
14
15# 4个cluster获取
16cluster <- read.csv('female_clustering.csv')
17female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
18table(female_clustering)

1.2 利用monocle V2构建对象

需要三样东西：表达矩阵、细胞信息、基因信息

1## 表达矩阵
2dim(female_count)
3## 细胞分群信息（包括6个stage和4个cluster）
4table(female_stages)
5table(female_clustering)
6## 基因信息
7gene_annotation <- as.data.frame(rownames(count_matrix))

开始构建

1library(monocle) 
2# 直接使用作者包装的函数，代码更简洁
3DE_female <- prepare_for_DE (
4  female_count, 
5  female_clustering, 
6  female_stages
7)

可以看到，包装好的函数只需要提供表达矩阵和细胞信息，其实它背后做的事情比较多，就像下面👇的代码

 1  # 1.表达矩阵
 2  expr_matrix <- as.matrix(count_matrix)
 3  # 2.细胞信息
 4  sample_ann <- data.frame(cells=names(count_matrix), 
 5                           stages=stage, 
 6                           cellType=cluster)
 7  rownames(sample_ann)<- names(count_matrix)
 8  # 3.基因信息
 9  gene_ann <- as.data.frame(rownames(count_matrix))
10  rownames(gene_ann)<- rownames(count_matrix)
11  colnames(gene_ann)<- "genes"
12
13  # 然后转换为AnnotatedDataFrame对象
14  pd <- new("AnnotatedDataFrame",
15            data=sample_ann)
16  fd <- new("AnnotatedDataFrame",
17            data=gene_ann)
18
19  # 最后构建CDS对象
20  sc_cds <- newCellDataSet(
21    expr_matrix, 
22    phenoData = pd,
23    featureData =fd,
24    expressionFamily = negbinomial.size(),
25    lowerDetectionLimit=0.5)
26
27  sc_cds <- detectGenes(sc_cds, min_expr = 5)
28  sc_cds <- sc_cds[fData(sc_cds)$num_cells_expressed > 10, ]
29  sc_cds <- estimateSizeFactors(sc_cds)
30  sc_cds <- estimateDispersions(sc_cds)

1.3 进行差异分析

其实也是包装了differentialGeneTest函数

 1start_time <- Sys.time()
 2female_DE_genes <- findDEgenes(
 3  DE_female, 
 4  qvalue=0.05
 5)
 6end_time <- Sys.time()
 7end_time - start_time
 8
 9# 1] "4435 significantly DE genes (FDR<0.05)."
10# [1] "3268 significantly DE genes (FDR<0.01)."
11# Time difference of 1.796054 mins

我们最后保留的也就是4435 significantly DE genes (FDR<0.05).

实际上，这个包装的函数就是下面这样，我们只需要提供创建好的monocle对象和q值即可：

 1function(HSMM=HSMM, qvalue=qvalue){
 2    diff_test_res <- differentialGeneTest(
 3        HSMM,
 4        fullModelFormulaStr="~cellType",
 5        cores = 3
 6    )
 7
 8    sig_genes_0.05 <- subset(diff_test_res, qval < 0.05)
 9    sig_genes_0.01 <- subset(diff_test_res, qval < 0.01)
10    # 帮我们统计了不同q值得到的差异基因数量
11    print(paste(nrow(sig_genes_0.05), " significantly DE genes (FDR<0.05).", sep=""))
12    print(paste(nrow(sig_genes_0.01), " significantly DE genes (FDR<0.01).", sep=""))
13
14    diff_test_res <- subset(diff_test_res, qval< qvalue)
15
16    return(diff_test_res)
17}

最后得到这么一个差异基因数据框：

上面得到的4435个差异基因，我们要知道这些基因在哪个cluster：

作者的做法是：先得到了每个差异基因在不同cluster的平均表达量，然后找平均表达量最大的那个cluster，就认为这个基因属于这个cluster

1# 有了差异基因以后，继续使用RPKM值
2load('../female_rpkm.Rdata')
3de_clusters <- get_up_reg_clusters(
4  females, 
5  female_clustering, 
6  female_DE_genes
7)

2 使用ROTS进行差异分析

2.1 准备表达矩阵和分群信息

 1rm(list = ls()) 
 2options(warn=-1) 
 3options(stringsAsFactors = F)
 4
 5# 6个发育时期获取
 6head(colnames(female_count))
 7female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1)
 8names(female_stages) <- colnames(female_count)
 9table(female_stages)
10# 4个cluster获取
11cluster <- read.csv('female_clustering.csv')
12female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
13table(female_clustering)

2.2 使用ROTS

下面会使用ROTS包（Reproducibility-optimized test statistic），对每个群和其他几个群的共同体进行比较

• ROTS可以使用RPKM值；它的速度可能会很慢

• 差异分析重点就在：表达矩阵和分组信息

配置最重要的分组信息

 1library(ROTS)
 2library(plyr)
 3# 使用RPKM矩阵
 4load('../female_rpkm.Rdata')
 5
 6# 当前分组是这样
 7> table(female_clustering)
 8female_clustering
 9 C1  C2  C3  C4 
10240  90 190  43 
11
12# 我们只需要用数字来表示，于是可以用substring(x, start, stop)获取
13groups<-female_clustering
14groups<-as.numeric(substring(as.character(groups),2,2))
15> table(groups)
16groups
17  1   2   3   4 
18240  90 190  43 
19
20# 如果要对第1群进行差异分析，那么就要把它和其他群比较(可将其他亚群定义为234)
21groups[groups!=1]<-234
22# 同理，对第2群进行差异分析
23groups[groups!=2]<-134

然后配置表达矩阵

1RPKM.full=females
2ROTS_input<-RPKM.full[rowMeans(RPKM.full)>=1,]
3ROTS_input<-as.matrix(log2(ROTS_input+1))

第1群

1start_time <- Sys.time()
2results_pop1 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)
3summary_pop1<-data.frame(summary(results_pop1, fdr=1))
4end_time <- Sys.time()
5(end_time - start_time)
6# Time difference of 18.40462 mins

然后第2群、第3群、第4群

 1groups[groups!=2]<-134
 2start_time <- Sys.time()
 3results_pop2 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)
 4summary_pop2<-data.frame(summary(results_pop2, fdr=1))
 5end_time <- Sys.time()
 6(end_time - start_time)
 7
 8# 第3群
 9groups[groups!=3]<-124
10start_time <- Sys.time()
11results_pop3 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)
12summary_pop3<-data.frame(summary(results_pop3, fdr=1))
13end_time <- Sys.time()
14(end_time - start_time)
15
16# 第4群
17groups[groups!=4]<-123
18start_time <- Sys.time()
19results_pop4 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)
20summary_pop3<-data.frame(summary(results_pop4, fdr=1))
21end_time <- Sys.time()
22(end_time - start_time)
23
24# 都得到以后，共同保存
25save(summary_pop1,summary_pop2,summary_pop3,summary_pop4,
26     file = 'ROTS_summary.Rdata')

“

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

”

---